Diet optimization problem

1. Introduction

Although in the modern world dieting is a hot and widespread topic due to the implications it has on health (obesity pandemic) and aesthetics (especially as social media constantly convey unrealistic expectations around body image), this blog post will reduce its discussion and analysis to the mathematical formulations that can arise from the challenge of choosing what and how much to eat given different criteria/objectives.

My first encounter with this problem was in Boyd's book (pg. 162), where the objective is to build a diet that contains enough nutrients to be a healthy one but is the cheapest possible. What does the "but" in the last sentence for? Because the cheapest diet without restrictions is to eat nothing 🤡.

For this study, the dataset used comes from this page, and the daily intake can be found in ./data/daily_intake.csv on the code repository.

You can take a look at the data on the following cell:

import pandas as pd

df = pd.read_csv('./data/stigler.csv')
print("df:")
display(df)

daily_intake = pd.read_csv('./data/daily_intake.csv')
print("daily_intake:")
display(daily_intake)

df:

	commodity	unit	price_cents	calories_k	protein_g	calcium_g	iron_mg	vitamin_a_kiu	vitamin_b1_mg	vitamin_b2_mg	niacin_mg	vitamin_c_mg
0	Wheat Flour (Enriched)	10 lb.	36.0	44.7	1411	2.0	365	0.0	55.4	33.3	441	0
1	Macaroni	1 lb.	14.1	11.6	418	0.7	54	0.0	3.2	1.9	68	0
2	Wheat Cereal (Enriched)	28 oz.	24.2	11.8	377	14.4	175	0.0	14.4	8.8	114	0
3	Corn Flakes	8 oz.	7.1	11.4	252	0.1	56	0.0	13.5	2.3	68	0
4	Corn Meal	1 lb.	4.6	36.0	897	1.7	99	30.9	17.4	7.9	106	0
...	...	...	...	...	...	...	...	...	...	...	...	...
72	Chocolate	8 oz.	16.2	8.0	77	1.3	39	0.0	0.9	3.4	14	0
73	Sugar	10 lb.	51.7	34.9	0	0.0	0	0.0	0.0	0.0	0	0
74	Corn Syrup	24 oz.	13.7	14.7	0	0.5	74	0.0	0.0	0.0	5	0
75	Molasses	18 oz.	13.6	9.0	0	10.3	244	0.0	1.9	7.5	146	0
76	Strawberry Preserves	1 lb.	20.5	6.4	11	0.4	7	0.2	0.2	0.4	3	0

77 rows × 12 columns

daily_intake:

	nutrient	daily_recommended_intake	unit
0	Calories	3.0	k Calories
1	Protein	70.0	g
2	Calcium	0.8	g
3	Iron	12.0	mg
4	Vitamin A	5.0	KIU
5	Vitamin B1	1.8	mg
6	Vitamin B2	2.7	mg
7	Niacin	18.0	mg
8	Vitamin C	75.0	mg

As data preprocessing steps we took df and:

• Split the "unit" column between "qtd" and the true "unit";
• Normalized the rows to contain the price and nutrients for 1 unit value of the true "unit" column;
• Applied accumulated inflation between 1939 and 2022 (this site was used). Although food price differences between 83 years are more complex than applying inflation, we used this simplification just to have an idea of what the price of the diets would be nowadays.

If you want to see the code go to the GitHub repository.

from dietopt.preprocess import apply_preprocess

df_processed = apply_preprocess(df)
print("df_processed:")
df_processed

df_processed:

	commodity	unit	USD	calories_k	protein_g	calcium_g	iron_mg	vitamin_a_kiu	vitamin_b1_mg	vitamin_b2_mg	niacin_mg	vitamin_c_mg
0	Wheat Flour (Enriched)	lb.	7.718292	4.470000	141.100000	0.200000	36.500000	0.0	5.540000	3.330000	44.100000	0.0
1	Macaroni	lb.	3.022998	11.600000	418.000000	0.700000	54.000000	0.0	3.200000	1.900000	68.000000	0.0
2	Wheat Cereal (Enriched)	oz.	5.188407	0.421429	13.464286	0.514286	6.250000	0.0	0.514286	0.314286	4.071429	0.0
3	Corn Flakes	oz.	1.522219	1.425000	31.500000	0.012500	7.000000	0.0	1.687500	0.287500	8.500000	0.0
4	Corn Meal	lb.	0.986226	36.000000	897.000000	1.700000	99.000000	30.9	17.400000	7.900000	106.000000	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...
72	Chocolate	oz.	3.473231	1.000000	9.625000	0.162500	4.875000	0.0	0.112500	0.425000	1.750000	0.0
73	Sugar	lb.	11.084325	3.490000	0.000000	0.000000	0.000000	0.0	0.000000	0.000000	0.000000	0.0
74	Corn Syrup	oz.	2.937239	0.612500	0.000000	0.020833	3.083333	0.0	0.000000	0.000000	0.208333	0.0
75	Molasses	oz.	2.915799	0.500000	0.000000	0.572222	13.555556	0.0	0.105556	0.416667	8.111111	0.0
76	Strawberry Preserves	lb.	4.395138	6.400000	11.000000	0.400000	7.000000	0.2	0.200000	0.400000	3.000000	0.0

77 rows × 12 columns

As now we have the data, the mathematical formulation may be clearer.

Our diet must contain $m$ nutrients with at least $b_1$, . . . , $b_m$ quantities, as you may see on the daily_intake dataframe we have $m=9$ and our vector $b$ will be the values in column daily_recommended_intake of daily_intake dataframe. To build the diet we choose nonnegative quantities $x_1$, . . . , $x_n$ of $n$ different foods, in our case $n=77$ which is the number of rows of our df_processed dataset. One unit quantity of food $j$ contains an amount $a_{ij}$ of nutrient $i$ and has a cost of $c_j$ (note that the values in df_processed are the transpose of the following matrix $A$). We want to determine the cheapest diet, $x_{opt}$, that satisfies the nutritional requirements, daily_recommended_intake. This problem can be formulated as the linear program (LP):

$$ \begin{equation} \tag{1} \begin{split} \text{minimize } \ & c^tx \\ \text{subject to } \ & Ax \succeq b \\ & x \succeq 0 \end{split} \end{equation} $$

Where if $a, b \in \R^n_+$, $a \succeq b$ means that $\forall i = 1, . . . , n \implies a_i \ge b_i$.

2. Finding the optimal diet

Now that we have the data and the mathematical formulation, we can model the problem and find the "optimal" diet given our criteria (the objective function). To do so, we will use CVXPY, a python package for solving convex problems. Although our problem is a linear one and CVXPY maybe feel like overkill, as linear problems are convex ones, this library is a more general tool than using a linear solver, meaning that we can keep using the same tool for other non-linear problems in the future.

The following code uses CVXPY to solve the problem of Eq. 1:

import cvxpy as cp

# Problem data
c = df_processed.USD.values
A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values

# Constructing the problem
x = cp.Variable(shape=len(c))
objective = cp.Minimize(c@x)
constraints = [
    A@x >= b,
    x >= 0,
]
problem = cp.Problem(objective=objective, constraints=constraints)

result = problem.solve()
print(f"{result = }")
print(f"{x.value = }")

result = 0.1536070376685601
x.value = array([-3.76573316e-14,  8.68348920e-14, -1.46738328e-13,  8.13624289e-14,
        5.34424634e-03,  2.95446061e-14,  3.06770764e-13,  4.37997030e-13,
        2.77853882e-13,  2.04462612e-13,  1.77817196e-13, -1.69845507e-14,
        6.44932975e-14,  3.49824397e-13,  9.76998354e-14, -7.88504668e-14,
        2.33506059e-14, -1.23905986e-14,  4.03105827e-14,  4.02896389e-14,
        1.10932580e-13, -1.13920983e-14, -4.78824654e-14,  2.04427572e-13,
       -3.59915089e-14, -2.24149401e-14, -1.19513449e-14,  2.23277626e-14,
        7.37335280e-14,  5.53354150e-14, -1.12987895e-15, -3.96683084e-14,
       -6.50090199e-15,  1.99263707e-14, -1.15558487e-14, -2.94876516e-15,
        4.50554631e-14, -1.87698517e-14, -4.19150510e-14, -2.37919029e-14,
        4.57123757e-13,  2.85788206e-13,  2.81256178e-14,  3.21319703e-14,
        2.43133822e-13,  1.13132451e-02,  4.68001567e-13,  1.77437607e-13,
        1.82422902e-13,  8.41262040e-13, -4.33849115e-14,  5.17574850e-03,
        2.39843674e-12,  2.17471565e-14, -4.93345153e-14,  2.21769476e-14,
       -1.03700236e-14,  1.33085121e-13,  6.30155777e-14,  1.06124977e-13,
        7.10903546e-14,  1.88576706e-13,  1.03854535e-13,  5.39432412e-14,
        1.95543157e-13, -5.32264168e-14,  3.87732822e-13,  2.89602726e-13,
        1.03066891e-01, -3.54310802e-14,  4.52763432e-14,  3.49870292e-14,
       -1.32270420e-13, -1.99661357e-13, -1.94840287e-13, -8.99133127e-14,
       -1.36720909e-13])

Translating the numbers, or the $x$ vector, of our problem, the optimal diet (weekly to make the numbers more readable) would be:

from dietopt.utils import compute_weekly_diet_df

df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
df_weekly_diet

	commodity	qtd_weekly	usd_weekly
0	Navy Beans, Dried	0.7215 lb.	0.91
1	Cabbage	0.0792 lb.	0.06
2	Spinach	0.0362 lb.	0.06
3	Corn Meal	0.0374 lb.	0.04
4		TOTAL	1.07

Ok ok, a diet that gives us enough nutrients and costs 1.07 USD/week seems pretty cheap, as was stated from the optimization problem. But, looking at the menu... will I just eat 4 different items? Yeah... this does not seem like a diverse diet that I'll be able to stick to. Can I model something like diversity on my diet?

x = cp.Variable(shape=len(c))
objective = cp.Minimize(c@x + 10*cp.norm(x,2)**2)
constraints = [x >= 0, A@x >= b]
problem = cp.Problem(objective=objective, constraints=constraints)
result = problem.solve()

df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
df_weekly_diet

	commodity	qtd_weekly	usd_weekly
0	Navy Beans, Dried	0.3731 lb.	0.47
1	Lima Beans, Dried	0.1857 lb.	0.35
2	Corn Meal	0.1935 lb.	0.19
3	Cabbage	0.0785 lb.	0.06
4	Spinach	0.0234 lb.	0.04
5	Peas, Dried	0.0169 lb.	0.03
6	Sweet Potatoes	0.0205 lb.	0.02
7		TOTAL	1.16

The little trick applied above was to add norm 2 regularization to the $x$ variable, the new optimization problem would be:

$$ \begin{equation} \tag{2} \begin{split} \text{minimize } \ & c^tx + 10\lVert x \rVert _2^2\\ \text{subject to } \ & Ax \succeq b \\ & x \succeq 0 \end{split} \end{equation} $$

But why the $l_2$-norm was applied? In a nutshell, if you derive the $\lVert x \rVert _2^2$ by some component of $x$, like $x_i$, you will find

$$ \frac{\partial}{\partial x_i} \lVert x \rVert _2^2 = 2x_i $$

Meaning that the bigger is $x_i$ biggest is the incentive (derivative) to decrease its value on the objective function. As $x_i$ gets close to zero the incentive decreases too, meaning, in practical terms, that we penalize big (relative) values and don't care too much about small ones. That would be something like: "I want a cheap diet but don't want to just eat one thing, as I want some variety, even if I can eat just a little".

Note: where does this random $10$ come from? It was chosen by try-and-error until a found a more diverse diet that didn't cost too much (find the balance between both objectives).

3. What was the problem again?

An interesting thing happened by the end of the previous section: we noted that our objective wasn't really aligned with what we wanted. This is an important conclusion that you might not see much value in, but we real-life we often see this misalignment between the mathematical formulation and the business need. When this misalignment happens, solving the optimization problem does not translate into solving the real-world problem.

So, initially, we wanted a diet that had all the necessary nutrients while being cheap as possible (see Eq. 1), but, what if, we were looking for getting lean, what would happen?

3.1 Trying to get shredded

This time we want to minimize our calories while ingesting enough of the other nutrients. The $\sigma$ vector represents how many calories each food contains.

$$ \begin{equation} \tag{3} \begin{split} \text{minimize } \ & \sigma^tx \\ \text{subject to } \ & Ax \succeq b_{\sigma} \\ & x \succeq 0 \end{split} \end{equation} $$

On equation 3, $b_{\sigma}$ is the modified version of our daily_intake where the minimum value for calories is $0$ (zero).

import numpy as np

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
b[0] = 0 # no need for calories constraint
calories = df_processed.calories_k.values

x = cp.Variable(shape=len(calories))
objective = cp.Minimize(calories@x)
constraints = [x >= 0, A@x >= b]
problem = cp.Problem(objective=objective, constraints=constraints)
result = problem.solve()

print(f"Calories: {result:.2f} kcal/day")

cost = df_processed.USD.values
df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
df_weekly_diet

Calories: 0.60 kcal/day

	commodity	qtd_weekly	usd_weekly
0	Salmon, Pink (can)	6.5861 oz.	18.36
1	Coffee	3.3034 lb.	15.86
2	Liver (Beef)	0.4617 lb.	2.65
3	Tea	0.4397 lb.	1.64
4	Celery	0.9029 stalk	1.41
5		TOTAL	39.92

It would be a more expensive diet but with a really small quantity of calories daily, only 600 cal.

Disclaimer: Please, remember this is a mathematical exercise, not a real diet, don't be dumb, and just eat salmon from now on.

3.2 Could we add "tastiness" to the mix?

Thinking about dieting we may realize that one of the big challenges of it is that people need to change their eating habits, meaning, they end up eating what they don't like. So, let's add a new column to our data set with the "tastiness" of the commodities, ranging from -5, "I really despise it", to +5, "If I could, I would eat this every day".

To avoid eating too much we also added an upper bound to the nutrients in the diet (which include calories, probably a big effect by my tastiness preferences 😁). The $\tau$ letter will represent the vector of the tastiness of each commodity.

$$ \begin{equation} \tag{4} \begin{split} \text{minimize } \ & -\tau^tx \\ \text{subject to } \ & Ax \succeq b \\ & Ax \preceq 1.2b \\ & x \succeq 0 \end{split} \end{equation} $$

PS: The negative signal on the objective function of Eq. 4 is to translate the tastiness maximization problem to a minimization one.

On the next cell, you may see the tastiness I see on the commodities available:

tastiness_df = pd.read_csv('./data/tastiness_by_romulo.csv')
tastiness_df

	tastiness	commoditie
0	0	Wheat Flour (Enriched)
1	3	Macaroni
2	2	Wheat Cereal (Enriched)
3	3	Corn Flakes
4	0	Corn Meal
...	...	...
72	5	Chocolate
73	0	Sugar
74	-2	Corn Syrup
75	3	Molasses
76	3	"Strawberry Preserves"

77 rows × 2 columns

Now, the code for solving the problem:

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
tastiness = tastiness_df['tastiness'].values

x = cp.Variable(shape=len(calories))
objective = cp.Minimize(-tastiness@x)
constraints = [
    x >= 0, 
    A@x >= b, 
    A@x <= 1.2*b
]
problem = cp.Problem(objective=objective, constraints=constraints)
result = problem.solve()
print(f"{result = }")

calories = df_processed.calories_k.values
print(f"Calories: {calories@x.value:.2f} k/day")

df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
df_weekly_diet

result = -12.753855640167883
Calories: 3.60 k/day
Cost: 71.15 USD/week

	commodity	qtd_weekly	usd_weekly
0	Chocolate	13.4992 oz.	46.89
1	Pineapple (can)	1.2057 No. 2 1/2	5.51
2	Evaporated Milk (can)	3.1825 oz.	4.57
3	Leg of Lamb	0.6486 lb.	3.84
4	Salmon, Pink (can)	1.3280 oz.	3.7
5	Coffee	0.6892 lb.	3.31
6	Butter	0.3695 lb.	2.44
7	Pork Chops	0.1338 lb.	0.88
8	Spinach	0.0112 lb.	0.02
9		TOTAL	71.16

Wooww, 71.16 USD/week is kind of expensive, couldn't we add the cost to the objective function so we could control it a bit?

The new problem formulation would be:

$$ \begin{equation} \tag{5} \begin{split} \text{minimize } \ & -\tau^tx + c^tx\\ \text{subject to } \ & Ax \succeq b \\ & Ax \preceq 1.2b \\ & x \succeq 0 \end{split} \end{equation} $$

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
c = df_processed.USD.values

x = cp.Variable(shape=len(calories))
objective = cp.Minimize(-tastiness@x + c@x )
constraints = [x >= 0, A@x >= b, A@x <= 1.2*b]
problem = cp.Problem(objective=objective, constraints=constraints)
result = problem.solve()
print(f"{result = }")

calories = df_processed.calories_k.values
print(f"Calories: {calories@x.value:.2f} k/day")

df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
df_weekly_diet

result = -3.88875776195683
Calories: 3.07 k/day
Cost: 44.21 USD/week

	commodity	qtd_weekly	usd_weekly
0	Chocolate	7.0156 oz.	24.37
1	Corn Flakes	7.0783 oz.	10.77
2	Evaporated Milk (can)	4.9459 oz.	7.1
3	Tea	0.3031 lb.	1.13
4	Corn (can)	0.2390 No. 2	0.53
5	Liver (Beef)	0.0367 lb.	0.21
6	Cabbage	0.0714 lb.	0.06
7	Spinach	0.0180 lb.	0.03
8		TOTAL	44.2

Ok ok, 44.21 USD/week seems more reasonable. As an effect of adding the cost to our objective function we got rid of Pineapple although I like it (+5 in tastiness) it seems a bit expensive compared to others. Because we now consider cost, Liver also got into my diet although I'm not such a big fan (-2 in tastiness).

Now we have an interesting problem: as we have two objectives in the same objective function, they might compete. Do they compete? We can test it by optimizing for each objective separately:

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
c = df_processed.USD.values

x = cp.Variable(shape=len(calories))
constraints = [x >= 0, A@x >= b, A@x <= 1.2*b]

objectives_dict = {
    "Cost": cp.Minimize(c@x),
    "Tastiness": cp.Minimize(-tastiness@x),
} 
for tag, objective_function in objectives_dict.items():
    print(f"Solving for: {tag}")
    
    problem = cp.Problem(objective=objective_function, constraints=constraints)
    result = problem.solve()
    print(f"{result = }")

    calories = df_processed.calories_k.values
    print(f"Calories: {calories@x.value:.2f} k/day")

    cost = df_processed.USD.values
    print(f"Cost: {cost@x.value*7:.2f} USD/week")
    
    df_weekly_diet = compute_weekly_diet_df(df=df_processed, qtds=x.value.tolist())
    display(df_weekly_diet)
    
    print('\n')

Solving for: Cost
result = 0.3957707708368031
Calories: 3.00 k/day
Cost: 2.77 USD/week

	commodity	qtd_weekly	usd_weekly
0	Liver (Beef)	0.1466 lb.	0.84
1	Milk	0.3343 qt.	0.79
2	White Bread (Enriched)	0.2973 lb.	0.5
3	Corn Meal	0.2758 lb.	0.27
4	Lard	0.0712 lb.	0.15
5	Onions	0.1580 lb.	0.12
6	Cabbage	0.0571 lb.	0.05
7	Peanut Butter	0.0112 lb.	0.04
8		TOTAL	2.76

Solving for: Tastiness
result = -12.753855640167883
Calories: 3.60 k/day
Cost: 71.15 USD/week

	commodity	qtd_weekly	usd_weekly
0	Chocolate	13.4992 oz.	46.89
1	Pineapple (can)	1.2057 No. 2 1/2	5.51
2	Evaporated Milk (can)	3.1825 oz.	4.57
3	Leg of Lamb	0.6486 lb.	3.84
4	Salmon, Pink (can)	1.3280 oz.	3.7
5	Coffee	0.6892 lb.	3.31
6	Butter	0.3695 lb.	2.44
7	Pork Chops	0.1338 lb.	0.88
8	Spinach	0.0112 lb.	0.02
9		TOTAL	71.16

As you can see, we have different results (different $x_{opt}$) depending on the objective. We can visualize it in a 2D plane where each axis is one objective:

import plotly.express as px

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
c = df_processed.USD.values

x = cp.Variable(shape=len(calories))
constraints = [x >= 0, A@x >= b, A@x <= 1.2*b]

objectives_dict = {
    "Cost": cp.Minimize(c@x),
    "Tastiness": cp.Minimize(-tastiness@x),
} 
x_list = []
y_list = []
for tag, objective_function in objectives_dict.items():
    problem = cp.Problem(objective=objective_function, constraints=constraints)
    result = problem.solve()

    cost = c@x.value
    tastiness_ = -tastiness@x.value

    x_list.append(cost)
    y_list.append(tastiness_)


fig = px.scatter(x=x_list, y=y_list)
size = 600
fig.update_layout(
    xaxis_title='Cost',
    yaxis_title='- Tastiness',
    width=size, height=size
)
fig.update_traces(marker=dict(size=20))
fig.show()

As you may see, we don't have a "perfect solution" where there the left-most point is the downwards-most point too.

To analyze the trade-off between the objectives we may give different weights for each objective and explore the curve that lies between the graphs of the previous plot. Let's do this!

4. Finding the Pareto curve

For finding the Pareto curve we will weigh the objectives using the variable $\mu$, as you may see in Eq. 6.

$$ \begin{equation} \tag{6} \begin{split} \text{minimize } \ & -(1-\mu)\tau^tx + \mu c^tx\\ \text{subject to } \ & Ax \succeq b \\ & Ax \preceq 1.2b \\ & x \succeq 0 \end{split} \end{equation} $$

The variable $\mu \in \R$ will assume many values between $[0,1]$ so we can see the trade-off between the two objectives.

from tqdm import tqdm

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
c = df_processed.USD.values

mu_list = [0] + np.logspace(-5, 0, num=1_000).tolist()

x = cp.Variable(shape=len(calories))
constraints = [x >= 0, A@x >= b, A@x <= 1.2*b]

cost_list = []
tastiness_list = []
for mu in tqdm(mu_list):
    objective_function = cp.Minimize(-(1 - mu)*tastiness@x + mu*c@x )
    problem = cp.Problem(objective=objective_function, constraints=constraints)
    result = problem.solve()

    cost_ = c@x.value
    tastiness_ = -tastiness@x.value

    cost_list.append(cost_)
    tastiness_list.append(tastiness_)

100%|██████████| 1001/1001 [00:05<00:00, 197.82it/s]

import plotly.graph_objects as go
import plotly.io as pio

pio.renderers.default = 'plotly_mimetype+notebook'

fig = go.Figure(go.Scatter(x=cost_list, y=tastiness_list, mode='lines'))
size = 600
fig.update_layout(
    title='Pareto curve',
    xaxis_title='Cost',
    yaxis_title='- Tastiness',
    width=size, height=size,
)
fig.update_xaxes(range=[0, 11])
fig.show()

So, what does this line means? For every point in this line, there is no other feasible point that simultaneously reduces cost and -tastiness at the same time, meaning, choosing a point in this line means you had made a trade-off between cost and tastiness.

As one outstanding professor once said during my undergraduate studies "Life is a small blanket: you cover your head and you end up uncovering your feet, you cover your feet and you end up uncovering your head". Meaning, on multi-objective optimization problems in real life the objectives always end up competing, there is no "solution to rule them all", and you will often need to choose a trade-off between objectives.

Does the Pareto curve only exist in 2D? No, but for visualizing it we usually stop in 3D. As an example, let's add the minimization of calories as a third objective:

$$ \begin{equation} \tag{7} \begin{split} \text{minimize } \ & \mu_0 c^tx - \mu_1\tau^tx + \mu_2\sigma^tx\\ \text{subject to } \ & Ax \succeq b \\ & Ax \preceq 1.2b \\ & x \succeq 0 \end{split} \end{equation} $$

In Eq. 7 the variables $\mu_0$, $\mu_1$, and $\mu_2$ will assume many values between $[0,1]$ so we can build the Pareto surface.

import itertools

A = df_processed.iloc[:, 3:].T.values
b = daily_intake.daily_recommended_intake.values.copy()
b_no_calories = b.copy()
b_no_calories[0] = 0
calories = df_processed.calories_k.values
c = df_processed.USD.values

n_points = 25
mu0 = [0] + np.logspace(-5, 0, num=n_points-1).tolist()
mu1 = [0] + np.logspace(-5, 0, num=n_points-1).tolist()
mu2 = [0] + np.logspace(-5, 0, num=n_points-1).tolist()
weight_combinations = itertools.product(mu0, mu1, mu2)

x = cp.Variable(shape=len(calories))
constraints = [x >= 0, A@x >= b_no_calories, A@x <= 1.2*b]

def optimize_diet(mu0, mu1, mu2):
    objective_function = cp.Minimize(mu0*c@x - mu1*tastiness@x + mu2*calories@x)
    problem = cp.Problem(objective=objective_function, constraints=constraints)
    problem.solve()

    cost_ = c@x.value
    tastiness_ = -tastiness@x.value
    calories_ = calories@x.value

    return cost_, tastiness_, calories_

from joblib import Parallel, delayed

results = Parallel(n_jobs=-1)(
    delayed(optimize_diet)(mu0, mu1, mu2) 
    for mu0, mu1, mu2 in tqdm(list(weight_combinations))
)

cost_list = [tuple_[0] for tuple_ in results]
tastiness_list = [tuple_[1] for tuple_ in results]
calories_list = [tuple_[2] for tuple_ in results]

100%|██████████| 15625/15625 [00:21<00:00, 739.86it/s] 

import scipy.interpolate as interp

mesh_n_points = 200
x, y = np.meshgrid(
    np.linspace(min(cost_list), max(cost_list), num=mesh_n_points), 
    np.linspace(min(tastiness_list), max(tastiness_list), num=mesh_n_points),
)

z = interp.griddata(
    points=(cost_list, tastiness_list),
    values=calories_list,
    xi=(x, y),
    method='linear'
)
fig = go.Figure(data=go.Surface(
    x=x,
    y=y,
    z=z,
    opacity=0.95
))
fig.update_layout(
    title="Pareto 3D curve",
    height=1000,
    scene = dict(
        xaxis_title = "Cost (x)",
        yaxis_title = "- Tastiness (y)",
        zaxis_title = "Calories (z)",
    )
)
fig.show()

5. Final thoughts

Although this diet optimization problem is well known and used for academic purposes, my objective with this blog post is to show the reader that defining the objective and translating it mathematically are not always easy tasks (heard of perverse instantion?), and, as most problems in real life, we will have competing objectives, meaning, choosing a trade-off between the multiple objectives is added to this already complex task.

That's it! I hope this was informative and made you think about properly setting objectives in your optimization problems. 💪